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We consider the scaling solutions of Smoluchowski's equation of irreversible aggregation, for a 
non gelling collision kernel. The scaling mass distribution /(s) diverges as s""^ when s — > 0. r 
is non trivial and could, until now, only be computed by numerical simulations. We develop here 
new general methods to obtain exact bounds and good approximations of r. For the specific kernel 
Ki{x,y) = {x^'° + y^'^'f, describing a mean-field model of particles moving in d dimensions 
and aggregating with conservation of "mass" s — {R is the particle radius), perturbative and 
nonperturbative expansions are derived. For a general kernel, we find exact inequalities for r and 
develop a variational approximation which is used to carry out the first systematic study of r(d, D) 
for Kf) . The agreement is excellent both with the expansions we derived and with existing numerical 
values. Finally, we discuss a possible application to 2d decaying turbulence. 
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INTRODUCTION 

Aggregation phenomena are widespread in Nature. 
They have such an impact on material sciences, chem- 
istry, astrophysics, that a large amount of literature has 
been devoted to them In such dynamical processes, 

particles or objects as different in geometry and size as 
colloidal particles, galaxies, small molecules, vortices in 
fluids, droplets, polymers, can merge to form a new en- 
tity when they come into close contact or interpenetrate, 
through diffusion (Brownian coagulation |^,^), ballistic 
motion (ballistic agglomeration mM) , exogenous growth 
(droplets growth and coalescence ) or droplet deposi- 
tion 

One is usually interested in the evolution of the statisti- 
cal distribution of the "mass" s, a quantity characteristic 
of each particle, that is conserved in the coalescence pro- 
cess: it can be either the actual mass, the volume, the 
area, the electric charge, or any other physical quantity, 
depending on the underlying physics. 

A great progress was achieved when it was proposed 
I p^ and observed both in real experiments and in numeri- 
cal simulations that the mass distribution N{s, t) exhibits 
scale invariance at large time: 



S{t) 



S{t) - 



(0.1) 



The divergence of the mass scale S{t) bears on the obliv- 
ion of initial conditions and physical cut-off or discrete- 
ness, as does the diverging correlation length of critical 
phenomena: universality arises in dynamics as well, with 
new universality classes. 

The exponents z and (3 are easily derived from conser- 
vation laws and physical arguments, but in many cases a 
polydispersity exponent t defined by f{x) ~ x~'^ when 
a; — > is observed, whose value is nontrivial though uni- 
versal. The prediction of r is still a challenge. 



Except for a few (usually ID) exactly solvable model 
[]l3| , ^ , analytical results are still lacking. The most 
popular, and the earliest, approach to these aggregation 
problems is Smoluchowski's equation a master equa- 
tion |13] for the one-body distribution N{s,t): 



N{si,t)N{s - si,t)K{si,s - si) dsi 



dN{s,t) _ 1 
dt ^ 2 



- N{s,t) N{si,t)K{s,si)dsi (0.2) 
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where the aggregation kernel K{x,y) is symmetric and 
is characteristic of the physics of the aggregation process 
on a more or less coarse-grained level. Such kinetic equa- 
tions are usually derived within a mean-field approxima- 
tion, where density fluctuations are ignored. Mean-fleld 
approximation is expected to be valid above an upper 
critical spatial dimension. This dimension is usually 2 for 
reaction-diffusion models, but van Dongen showed that 
it can depend on the kernel JlGf . Including some proper 
approximation of the density-density correlations in the 
kernel may improve Smoluchowski's approach p7[ |. 

Mean-field as it may be, Smoluchowski's equation is 
still highly nontrivial. No exact solution is available, ex- 
cept in a very few specific cases (see below) , and extract- 
ing the nontrivial exponent r for a specific system from 
the proper kinetic equation is not an easy task. The 
problem was clarified by van Dongen and Ernst [ p^ who 
classified the kernels according to their homogeneity and 
asymptotic behavior: 



K{bx,by) = b^K{x,y) 
K{x,y) - x^^y" [y > x) 



(0.3) 
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For a given physical system, the homogeneity A is easily 
determined using scaling arguments. We consider only 
nongelling systems with A < 1 [Q. For /i > 0, the expo- 
nent T is trivial and found to be t = 1 + A, whereas for 
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/i = 0, T depends on the wh ole solution / o f th e scaling 
equation derived from Eq. (0^) (see Eq. (1^) below), 
/i < does not lead to any power law behavior but rather 
to a bell-shaped scaling function / [|l8| . 

In the following, we shall focus on the fj, — case for 
which the exponent r has so far only been determined 
numerically by direct simulation of Smoluchowski's equa- 
tion 21|, not an easy task by time series [22[, 
and of course by direct simulation of the physical sys- 
tem supposed to be described by the considered Smolu- 
chowski's equation |l^,|22l. In the latter case, di- 
rect comparison with mean-field results is in principle 
rather delicate. These methods are quite heavy, which 
explains that very few values of r are known pl| , p2[ , 
most of them concerning a specific kernel, Kf){x,y) = 
{x^^^ + y^/^)'^ {0 < d < D), which appears in various 
physical applications . 

Considering the ubiquity and the importance of the 
/J, = case leading to nontrivial polydispersity exponents, 
analytical results as well as more effective numerical 
methods, making it possible to carry out extensive stud- 
ies, are certainly needed to use Smoluchowski's approach 
in a predictive way. The purpose of this article is to pro- 
vide both and use them to perform the first complete 



study of r(d, D) for the kernel Kf^ = {x^^^ 
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These analytical methods consist of exact bounds, per- 
turbative and nonperturbative expansions around ex- 
actly solvable limits, while we introduce a variational 
scheme, leading to excellent approximations of r at ex- 
tremely low computational cost, without directly solv- 
ing Smoluchowski's equation. We end the paper with a 
practical application of our results in the field of two- 
dimensional turbulence. 

In section |, we present a mean-field model of ag- 
gregation of Z?-dimensional spheres diffusing in a d- 
dimensional space and coalescing with conservation of 
their volume, for which we derive a Smoluchowski's equa- 
tion with the kernel Kf^ = (x^^^ + y^/^^. Under the 
scaling hypothesis, we write down the equation for the 
scaling function, determine the exponents z and (3, and 
derive an integral equation for r as well as a series of in- 
tegral equations for the moments of the scaling function 
/. This section introduces hardly any new result and is 
intended merely to clarify notations, to present the state 
of the art, and to make a few useful remarks. 

Section |l| and III present new analytical results for the 
previously introduced kernel Kf^ . Section |l| describes a 
method to obtain exact bounds for any kernel, based on 
integral equalities established in section ||. 

Section [II deals with expansions of r around its value 



for exactly solvable kernels. Starting from the remark 
that Kf) reduces to the constant kernel in both d ^ 
and 13 — > oo limits, for which an explicit exponential so- 



lution is known, we find some perturbative expansions in 
both limits. In the large D limit with d/D = A fixed, the 
kernel reduces to 2'^{xy)^ and we show that t ^ 1 + A, 
the first correction being exponentially small at large d, 
and thus nonperturbative. 

In section IV, we present a variational approximation 
based on integral equations for the moments of /, and 
valid for any homogeneous kernel. This method repro- 
duces some known exact results, and is used to compute 
T for a wide range of d and D, the results being summa- 
rized on Fig. ^ The approximation is compared to the 
few existing numerical results l21^,p3l as well as with an- 



alytical expansions derived in section [II, with excellent 



agreement and very low computational cost. 

Section presents a possible application in the field 
of two-dimensional turbulence. We consider a model 
of diffusing and merging coherent vortices, and Smolu- 
chowski's equation leads to non Batchelor energy spectra 
with exponents in qualitative agreement with direct sim- 
ulations fomid in the literature i3G,Pli. 



I. MODEL AND SCALING 

Consider hyperspherical particles in a d-dimensional 
box, of polydisperse radii R with distribution F(R,t), 
evolving the following way: at time t we choose the posi- 
tions of their centers with uniform probability in d-space. 
Then each pair of overlapping spheres of radii Ri and i?2 
merges to form a new sphere of radius. 



(1.1) 



where D is a parameter with D > d. D can be the ac- 
tual dimension of the spheres, as for instance in the case 
of _D = 3 spheres deposited on a d = 2 plane |l|] . Once 
each coalescence has been resolved, we have reached time 
t + St. 



A. Derivation of Smoluchowski's equation 

The conserved variable is s = R^ , and is continu- 
ous. We shall call sq the physical lower cut-off, that 
is the charge of the smallest sphere in the initial con- 
dition. Since the radius of a surviving sphere can only 
increase through coalescence, N{s,t) = for s < sq and 
for any time t > 0. Smoluchowski's equation consists just 
in a balance of collisions. The number of collisions be- 
tween two spheres of radius and S2^^ randomly and 
independently deposited in the d dimensional medium 
being N{si,t)N{s2,t)nd{s\^° + sl^^^ where Via is the 
d-dimensional total solid angle. We obtain the equation, 



N{s, t + St)- 7V(s, t)^nd\- / N{si,t)N{s - sut)K%{si,s - si) dsi 
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-Nis,t) 
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N{si,t)Kf){s,si)dsi 



(1.2) 



with K%{x,y) = {x^/'^ +y'^/°Y . We can get rid of the 
multiphcative constant, by properly choosing the time 
unit 5t and by replacing the finite difference in time by a 
partial derivative to exactly obtain Eq. (0.2). We notice 
that the only approximation used to derive the equation 
is to neglect multiple collisions, for the system is intrin- 
sically mean-field. 

The kernel K%{x,y) = (x^/-° + y'^/^Y has been in- 
troduced in many contexts from molecular coagulation 
I prt to cosmology |^,^ for specific values of d and 
£), and is one of the most studied in the literature 
|0,||j2l]-|2|,|2|-||| although very few analytical results 
are known. This kernel has A = and /i = 0. Ex- 
act solutions are available in the case d = or D = oo 
(constant kernel) Q , and d = D = 1 . 



B. Scaling 

Now, we introduce the scaling form of N{s,t). We 
first write the conservation law. The total mass in the 
system is f^^°° sN{s,t)ds ^ S{t)^~^ J^°° xf{x)dx and 
is conserved which implies (3 = 2, implicitly assuming 
that the integral of xf{x) converges, i.e., in terms of 
the small x divergence of /, that r < 2, which will be 
shown below. We consider the total number of parti- 
cles in the system n{t) = J^°° N{s, t) ds. It behaves 
at large time like S{tY^^ ^'^Tg,^^ f{x)dx. If r < 1, 



r 



n{t) ^ S{ty ^ °^ f{x)dx whereas if r > 1, n{t) oc 



J 



S{tY-^. If T = 1, the integral diverges hkc ln(S'(t)), 
hence n{t) S{tY-P\nS{t). 

As promised, we are now able to show that r < 2. if 
T > 2, the total charge in the system is proportional to 
S{ty~^, enforcing /3 = r. As a consequence, n{t) would 
have a non zero limit which is impossible. To summarize 
these results, we have, with n(t) cx , 



z, 

z{2- 



if T < 1 

r), if T > 1 



(1.3) 
(1.4) 



The derivation of the scaling equation is rigorously de- 
scribed in jl^, where it is shown that S{t) ~ wt^, w 
being some positive constant characteristic of the time 
dependent equation. Plugging the scaling form of the 
distribution into Smoluchowski's equation, and match- 
ing the large t behavior of both sides of the equation, 
yields z = D/{d — D) and the equation for the scaling 
function. 



r+oa 

w [sfis) + 2f{s)] = fis) / f{s^)K%{s,, s)dsi 
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/(si)/(s - si)K'^isi,s - si) dsi 



(1.5) 







If r > 1 each term of the RHS of Eq. (L5) is separately 
divergent and they should be properly grouped, for in- 
stance. 



-t-oc 



w [sf'is) + 2f{s)] = f{s) I f{si)Ki{si, s)dsi 

Js/2 

f{si) [f{s - si)kUsi,s - Si) - f{s)Ki{sus)] dsu 



(1.6) 
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another way of taking care of these divergences is to 
be found in 

As we are only interested in the exponent affecting the 
small s behavior of /, we shall set w to unity by changing 
/tow/. If /(s) is a solution of Eq. (^), then &i+^/(fos) 
is also a solution. The value of b is often fixed by impos- 
ing Jxf{x)dx = 1, but we will make a different choice 
for reasons that will become clear later. 

A careful study of the large s behavior of / shows 

Coo S s~'^e~^'' , with 

^dx Hi. We choose 
the solution corresponding to S = 1, which fixes &, and 
leads to a nontrivial value for Jxf{x) dx. This asymp- 
totic behavior is not valid for A = 1 {d = D). 

For d = or £) = oo, Eq. (pTs]) reduces to the constant 



that if A < 1 (d < D), f{s) R 
--'-J^^'kUxA-x)x-'{1-x) 



kernel equation with exact solution fo{x) = 2e~'^ and 
/oo(s) = 2^^'^e^'' (note that the large s asymptotics be- 
come the exact solution for all s in these cases) . For d = 1 
and _D = 1, an exact analytic solution is also known for 
the time dependent equation, the scaling function being 



oc s 



-3/2, 



pa], with 2 = oo and S{t) cx e* 



Now, for given d and D, and plugging the expected 
small s behavior f{s) ^ s^'^ into Eq. (1.5), one first gets 
that r< 1-|-A = 1-|- d/D . Then, matching the behavior 
of both sides of Eq. 1^ |l|,|l9|, one finds, 



T = 2- 



x^dx. 



(1.7) 



If a > r — 1 we obtain by multiplying Eq. (1.5) by x° 
and integrating [p8|,p7[, 
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2(1 - a) / x°'f{x)dx 



f{x)J{y)Ki{x, y) K + y"- (a; + y)"] dxdy. 



(1.8) 
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C. Existing analytical and numerical results 

Most existing analytical results for /i = kernels are to 
be found in the beautiful series of papers by van Dongen 
and Ernst |T6|Jl^ , p^ , ^ . Apart from results mentioned 
earlier, they determined the small x subleading behavior 
of the scaling function, and they found some inequalities 
for T in the cases d = 1, and D = 1. In 1984, Leyvraz 
psf proposed the analytical result r = 1 + 1/2D for the 
kernel Kf^ with d = 1, but in 1985, using exact inequal- 
ities, van Dongen and Ernst showed that this result was 
erroneous and explained why it was so |p7| . The argu- 
ment of Leyvraz leading to this result is perfectly valid 
for class I kernels with ^ > for which it predicts the 
correct exponent, but it breaks down for ii — kernels. 
We mention this fact for some references to the wrong 
result T = 1 + 1/2D can still be found in some recent 
articles. 

We now review various kinds of numerical studies con- 
cerning the polydispersity exponent r. These studies deal 
with the kernel Kfy. 

Kang et al. simulated a model of particle diffusion 
and coalescence (PCM) that can be shown to be exactly 
equivalent to Smoluchowski's equation. They also nu- 
merically directly computed the solution of the equation 
itself. Their results concern the d = 1 case. They sur- 
prisingly found values of t in contradiction with the ex- 
act bound T > 1 (see section |^ (for D = 4, they found 
T — 0.63). By comparison between their two methods 
of computation, they concluded that in both cases they 
observed a pseudo-asymptotic state, with wrong expo- 
nents but apparent scaling, and that the actual asymp- 
totic scaling regime appeared at times too large to be 
seen by their simulations. This illustrates the drawback 
of considering the direct time evolution of the system: 
the actual asymptotic regime may not be reached within 
the accessible numerical simulation time scale. 

Krivitsky ||2l[] numerically solved Smoluchowski's 
equation for the time dependent distribution for the ker- 
nel Kf), for Z? = 1, c? < 1, for which he determined 10 
values for r (see Fig. |^). Comparison with analytical re- 
sults obtained by analysis of the scaling equation (infinite 
time limit) in the present article will assess the fact that 
in this case the asymptotic regime was actually reached 
by Krivitsky's solution. These numerical results will be 
found to be in excellent agreement with our variational 
method of section pv|. 

Song and Poland ||22(|, computed the large time evo- 
lution of the number of clusters n(t) oc t^^ , and as 
z' = z(2 — t) when t > 1, and z' = z, when t < 1, 
we can extract r from their data (for which t > 1). 
Their method consists in solving the equation for n(t) 



as a power series in time t, and to extract the expo- 
nent z' by manipulations of this series. They treated 
only the cases d = 1, D = 2 and d — 2,D — 3. In 
the case d — 1,D ~ 2, they present two different re- 
sults in the text (?). They first consider K2 and find 
1/z' = 0.57 ± 0.01, then they extend their method two 
K^^^ and in the case d = 2, which is exactly the same as 
previously, they find 1/z' = 0.588 (they do not give any 
error estimate in this case). In the following, we shall see 
that we believe the first result to be closer to the exact 
one. In the next section, we shall see that their result 
in d = 2, D = 3 strongly violates exact inequalities, and 
thus is wrong. 

The conclusion of this section, is that no complete 
study of the value of r had been performed until now 
because of a lack of appropriate numerical tools. More 
precise analytical results would also certainly be welcome 
to guide numerical works. We see that simulating or solv- 
ing for the time evolution of the distribution function may 
not enable to reach the asymptotic scaling regime, and a 
guideline of the present work will be to directly rely on 
the scaling equation corresponding to the infinite time 
asymptotic state itself. 



II. EXACT BOUNDS 



In the next three sections, our workhorses will be both 



Eq. (|17D and (|1J). 

We first show that r > 1, for d > 1. Suppose r < 1 
and consider Eq. ( |l.8| ) with a — 0, 

r+oc p p+ao 

2 f{x)dx= fix)f{y)KUx,y) dxdy. (2.1) 

Jo J Jo 

For d > 1, we have {xi +y^Y > x'i -f y'^', which leads 
to Jf{x)dx > Jf{x)dx J f{x)xodx (in the bulk of the 
text, all integrals should be understood from to 00). 
Comparing with Eq. ( |1.7| ), this leads to 1 > 2 — t o r 
T > 1, which is contradictory. Notice that Eq. ( |l.8| ) 
with a = 2 for d = 1 and D = 1 leads to / x'^f{x)dx — 
2{Jx'^f{x)dx){Jxf{x)dx), and we recover the exact re- 
sult T — 2 — Jxf{x) dx — 3/2 |2^ in a very simple way. 
These results were already obtained by van Dongen and 
Ernst 1 27 1^ , who were able to find in the case D = 1 the 
exact inequality, 2d < t < 2-2i-^(l-d)/(2-2''), which 
shows that r = 2d-|-0(d^) when d 0. This interesting 
result will be generalized to any D in next section and 
the O(d^) term will be computed m D = 1. They also 
found weaker inequalities in d = 1, but no result was 
obtained for general d and D. 
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In order to deal with a general ^ = kernel K, we in- 
troduce an extremely simple meth od t o get lower and up- 
per bounds for r. We rely on Eq. ( |l.8| ) valid for a > r — 1 . 
Combining Eq. (L7) and (lOI) , we get: 



9{x,y)dxdy 



(2.2) 



where A(w) = {l + u" - {l + u)'')K{x,y)/{u°' +u^) satis- 
fies A('u) = A{l/u) and qix,y) = {x"y^+x^y°')f{x)f{y). 
The ratio in Eq. ( |2.2[) can then be interpreted as the in- 
verse of a kind of average of A{x/y) with the weight 
g{x, y). For a given a < A, we numerically determine the 
maximum Mg and minimum nia of the function A(u). 
Using Eq. (]2.2D, this gives 



(2.3) 



2 - (1 - a)/TOa < r < 2 - (1 - a)/Ma 



We then choose the best values of a < A compatible 
with a > T — 1 leading to the tightest bounds. More pre- 
cisely, we proceed the following way: we start with a = A 
(as r < 1 -f A), from which we obtain some upper and 
lower bound r,„ and tm- If tm < 1 -|- A, Eq. (2.3) holds 
for Tm — 1 < a < a, and we can compute new bounds 
for each a in this interval, and find the tightest bounds. 
The upper bound obtained for a — \ cannot be improved 
since A{Q) = 1 for a < A, hence 2 - (1 - a) /Ma > 1 + a, 
but in many cases we can find a better lower bound. 

For _fC^, a superficial plot of the function A{u) may 
lead to the incorrect conclusion that its minimum is al- 
ways obtained at m = with A{Q) — 1. In fact a more 
careful study of A shows that for certain values of a, the 
actual minimum is at u > but very close to 0. For 
M ^ 0, A{u) 1 + du^/^ - m'*/^-", and we see that 
if a > (d — ^)/D, there is a local minimum for Um > 
with A{um) < 1. For d > 1, and a = (d — 1)/D + e, we 
get Um ^ exp(— ln(c?)/e), which vanishes exponentially 
when £ — > (d > 1). Indeed, even when a is not so close 
to {d — 1)/D, Urn may be very small. For instance, for 
d = 2, D = 3, and a = 0.58598 > {d - l)/D = 0.333..., 
we find that Um = 1-365 x lO"'', and A{um) = 0.7322, 
which leads to a nontrivial lower bound of 1.4349 for t. 

Actually, it is easily seen that the inequalities obtained 
by van Dongen and Ernst (in the case d = 1 ox D = 1) 
correspond toa — d/ D . In fact even in this case, and 
ma are nontrivial, and they used some explicit bounds of 
Ma and m^, which do not lead to the tightest bounds 
for T. 

Thus, our method consists in computing the actual 
value of nia and Af^, and varying a to optimize these 
bounds, which allows us to greatly improve van Dongen 
and Ernst's explicit inequalities for _D = 1 or d = 1, 
and to obtain new exact bounds for d > \. For in- 
stance, for the physically interesting cases (see below) 
(d = 1, D = 2), [d = I, D = A) and {d = 2, D ^ 4) 
we respectively found 1.084 < r < 1.147, 1 < r < 1.075 
(compared to 1 < r < 1.28 and 1 < r < 1.109 in (27)) 
and 1.25 < r < 1.5. 



For d = 2,D = 3, we find 1.4349 < r < 1.585, which 
just discards the value r = 1.244 found by Song and 
Poland and strongly questions the validity of their 
approach. The exact bounds we obtained in d = 1, 1? = 2 
are violated by their alternative value 1.1 50 fo r r but not 
by their first result 1.123 (see subsection I C). 

It is useful to note that for any D, with a — d/D, 
A{u) — » 1/2 when d ^ 0, which entails that r — > (from 



Eq. (gj)) in this Hmit. 

To conclude with this topic of inequalities, let us con- 
sider Eq. (2.3) with a — d/D. In this case, when 

D —>■ CXD, 



1 



A{u) = -{1 + u-oy 



1 + UO - [1 + u)o 



d-1 



<u<l 
= 



(2.4) 



hence nia — > 1/2 a nd M a — > 2"^ ^. Therefore, the upper 
bound for T in Eq. (U) tends to 2- 2^-''. This is strictly 
less than 1 for d < 1, which means that for any d < 1, 
there exists a finite critical Dc{d), such that r < 1 for 
any D > Dc- This result will be used in section [II. 



III. PERTURBATIVE AND 
NONPERTURBATIVE EXPANSIONS 

In this section we use the exactly solvable limits d = 
and D = oo as a basis for a perturbative expansion. We 
also consider the case d —>■ go, keeping d/D = X constant, 
for which we find a nonperturbative expansion. 

We saw that limd^o t = 0. What about the D ^ oo 
limit of T ? In fact, although strictly at D = oo, 
T is equal to 0, as f{x) — 2^~'^e~^, we will see that 
Too = hm£)_+oo T > 0. This result was already noticed by 
van Dongen and Ernst in d = 1 Since t < 1 + d/D 
we get that. 

Too < 1 (3.1) 

What can we learn from equation ( |l.7| ) in the large D 
limit ? We see that the limit for r is 



Too 2 - / foo{x) dx = 2 - 2 
/o 



(3.2) 



provided that: 



lim 



+ 00 



ifDix)^fooix))xodx = (3.3) 



For d < 1, this result is consistent, since, from the last 
remark of section ||, we get Too<2 — 2^^'*<1. 

However, for d > 1 we know that t > 1, hence Tqo = 1, 
which means that for d > 1, 



lim 



+ 00 



ifoix) - foo{x))x* dx=l- 2^-'' > 



(3.4) 
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while in d = 1, (3.3) is true. 

Now that we know the large D limit of r (tqc = 1 
for d > 1 and Tqo = 2 — 2^^'' for d < 1), as well as its 
small d limit (r — > 0), let us compute the corresponding 
asymptotic corrections. 



A. Small d expansion 



First, consider the limit d ^ 0. We expand / in series 
in d: fix) = /o(x) + + 0{d^), fo{x) ^ e-\ A 

systematic way of expanding r would be to write down a 
linear (self-consistent) differ ent ial equation for /i to solve 
it and plug the result into (1.7). 

However, as far as the first order is concerned we can 
get it without solvin g for /i. By developing the integral 
expression of t, Eq. (1.7), we get. 



f{x)x o dx 



d r°° f 
15 I fo{x)lnxdx - d I 

+0{d^). 



+ OC 



fi{x)dx 



(3.5) 



Then we expand both sides of Eq. (21) to get an equa- 
tion for J fi{x)dx, 



-i-oo 



fi{x)dx 



+ 00 



foix)fo{y)Hx^ +y^)dxdy 



+ 00 



fn{x)dx I fi{x)dx 



(3.6) 



hence Jfi{x)dx = ^//e ^ ^ln(a;D + yo^dxdy. After 
eliminating J fi{x)dx, we get: 



r = 2dJD + 0{d^ 



Jd 



In 1 



1-u 



du 



(3.7) 



Let us mention that this result can be systematically gen- 
eralized to the case of any homogeneous kernel of the 
form: {g{x,y)Y, leading to, t = 2d lng(l, i^)<iw -|- 

o{<P). 

Although it may seem a bit tedious, it is interesting 
to recover this result in another way, as it shows that 
the small x behavior of /i is consistent with the d ^ Q 
expansion of the power law x""^ = 1 — 2dj£i Inx + 0{d^). 
Let us write down the linear equation for /i , 



xf[{x) + 2e-- / /i(y)e«dy = 2e-^ / h{y)dy + 4e 



+ CX3 



'nn{y^''^ +x^''^)dy 



2e~- / Hy^'^ + {x~yf'^)dy. 



(3.8) 



With u = e^/i we get the following equation: 



x{u' — u) + 2 / u{y)dy 



-t-oo 



u{y)e ^dy + 4 



-t-oo 



'ln(?/ 



1/^ I ^i/D 



)dy 



2xJd — —{xlnx — x), 



(3.9) 



which implies, after taking the derivative of Eq. (3.9), 



xu" + (1 — x)u' + " ^ ^ 



X 



1/D-l 



y 



1/D , ^1/D 



dy 



D 



\nx- 2Jr 



(3.10) 



the solution u of (|3.10|) involves two integration con- 
stants, one being fixed by the fact that /i should go to 
zero at large a;, the other, co, by writing the compati- 
bility with Eq. (3.9), which can be done by taking the 
X — » limit the latter equation. From the expression of 
the solution (appendix HI), or directly from Eq. (3.10), 
it is easily seen that u has the asymptotic expansion for 
a; ^ 0: 



with bo = cq — 2/D. 

We know that f{x) ^ cx~'^ when x ^ 0. When d ^ 0, 
c — ^ 2 and t — dri + 0{d^), hence up to order d we ex- 
pect. 



f(x) ~ 2ti Inx 
so that we interpret &o as — 2ti, 



T = — rf 

2 



0(rf2 



The .X ^ limit of (3.9) is : 



+ 00 



fi{x)dx 



4 
D 



lux dx 



(3.12) 



(3.13) 



(3.14) 



u{x) = bolii-x + 0(1) 



(3.11) The integration of Eq. ( [3.1C ) between and -l-oo yields. 
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+ 00 



-bo+ fi{x)dx ^ 2Jd - 



D 



e ^ In a; da; 



T = 2d+{—-A)(f + 0{d?) 
o 



(3.16) 



(3.15) 

The combination of Eqs. (3.14) and ( 3.15| ) yields 6oi 



In section |^ (see Fig. |^), we shall see that this result is 
in excellent agreement with both Krivitsky's results and 



which, substituted into Eq. (flsl), eventually leads to ^ ^^w method of approximation that we shall introduce 
the same result for r as previous ly o btained through the ™ section ^ . 
expansion of Eq. (^.l|) and Eq. (1.7). 

For Z? = 1, we get t — 2d + 0{d^), in good agree- 
ment with direct numerical integration of Smoluchowski's 
equation performed by Krivitsky [ pT[ and shown on Fig. 
|l| (see below). This result for D ^ 1 also coincides up to 
order 0{d) with the inequalities for t that we obtained 
above, as noticed in section ||. This is not the case for 
other values of D. 

The order 0{d^) requires the computation of /i. How- tioned in section M, r < 1 for any D a bove a finite critical 
ever, in the special case Z) = 1 it is possible to obtain Dc{d). As a consequence, Eq. (L8) can be written for 
explicitly the 0{d^) term by expanding Eq. (1.8) for ai\y D > Dc{d). Therefore, we can develop this equation 



B. Large D expansion 

Now, we perform an expansion in powers of 1/D for 
d < 1, expanding f{x) = foo{x) -f -^/i(a;) -I- -^hix) + 
Oil/D"^). 

Perturbative expansion m d < 1 - In d < 1, as men- 



a — d/D (see appendix ||). We obtain 



J 



for large D in powers of and we find at first order. 



+ CXD 



fi{x)dx = d2 



d-2 



/oo(a;)/oo(y)(lna; + Iny) dxdy + 2" 
I 



+ 00 



+ CXD 



foo{x)dx / fi{x)dx, 



(3.17) 



hence 



+ 00 



fi{x) dx 



+ C30 



foc{x) In (a;) da; 



We conclude, using Eq. (3.18), that the first order cor- 
rection to Too is zero. 

The same method also gives access to the next term: 



(3.18) 



where 7 is Euler's constant, while from Eq. (1.7), 
T = Too - — y fooix) lii{x) dx 



2D^ 



+00 



/oo(a;)(lna;) dx 



d 



+00 2 
fi{x) hi{x) dx — 



+ 00 



/2(a;)dx, (3.20) 



DJo 



+ CX2 



f,{x)dx + 0{^). 



(3.19) 



while, 



J 



+ 00 



f2{x)dx = - 



+00 

foo{x)foo{y)^ [(d-t- l)((lna;)2 + (Iny)^) -f 2(d - 1) ln(x) ln(y)] dxdy 



h {x)h (y) (In x + h\y) dxdy 



-\-OQ 



+ 00 



h{x)dx] +2 / /2(a;)da; 



Using the known value of J /i and our favourite integral table, we get: 



f2{x)dx 



-\-oo 



+ 00 



foo{x){\nx) dx + d fi{x)\n{x)dx + 



(3.21) 



(3.22) 



pansion for t without solving for /i and /2 themselves, 
although this can also be achieved this way . Note that 
_ 9 _ 91-d I TT 2 d(l — d) _}_\ (n no^ iu thc limit of large D and small d, Eq. (3.7) and ( 3.23| ) 
^" + ^--^ K^-^-^) coincide up to order 0(d/i?2). 

Perturbative estimate for d > 1 - In the case d > 1, we 



(7 being Euler's constant), which leads to: 
Once again we were able to obtain a highly nontrivial ex- 
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have shown that r > 1 and since r < 1 + d/D, we see 
that T — » 1 for D — >■ oo and finite d > 1. As /i is non 
integrable, Eq. (1.8) cannot be used with a — 0, and the 
previous perturbation breaks down. 

Nevertheless we can try to obtain an estimate of r 
in the following way: we make th e an satz / foo + 
c/s^+^e-". We plug it into Eq. ([r|) and Eq. ([ij) 
for a — d/D, and after some algebra (see appendix 
[V ) we see that for consistency e must be of order 1/D 
and that c (1 — 2^~'^){d/ D — e), and eventually that 
e = k/D + 0{1/D^) where k is the solution of the non- 
linear equation: 



1 



{l + v^)'^dv 



(3.24) 



This equation always has a solution consistent with the 
exact bound 1 < t < 1 + d/ D. For instance in the case 
d = 2, D = 4 we obtain r « 1.462. Though it is stiU of 
order l/D, the obtained perturbative estimate depends 
on the choice of a. a ~ d/D seems however to be the 
most natural choice. 

In d = 1, c vanishes and we do not learn much. All 
terms of the d < 1 series for t in powers of 1 /£> van- 
ish for d 1, as can be seen in Eq. ( 3.23 ) for the 
two leading ones. The reason is the following: the per- 
turbation is derived from Eq. ( ^.l[ ) under the assump- 
tion that T < 1. In d = 1, such an assumption yields 
2 / f{x)dx = 2{J x^/"f{x)dx)ij f{x)dx) hence r = 1. 
Consequently the perturbative value of r tends to 1 when 
d ^ 1~. As will be illustrated below by numerical re- 
sults, for a given d > 1 there is a critical D = Dc{d) 
above which r < 1, and Dc{d) tends to infinity when 
d — > 1~, entailing the vanishing of the perturbation va- 
lidity domain in D. Thus, the correction to r = 1 for 
large D may be nonperturbative in d = 1. 

If we now take the d — > oo limit in Eq. (3.24), we 
obtain, r ~ 1 + A — (A = d/D), a nonperturbative 

behavior in d which is to be related to the results below, 
obtained for d — >• oo, D — » oo, keeping A constant. 



Kix, y) = 2'^(xy)^ [l + 2-^0(l/d2)] (3.26) 

The rescaled function / = 2*^/ is the solution of the 
scaling form of Smoluchowski's equation with the ker- 
nel K = 2^'^K{x,y), which is equal to {xy)^ at every 
order in 1/d — \/{XD). In fact, we can estimate this 
correction by assuming that for finite d and D, 



f{s) - cx/s 



(3.27) 



for s 0. Plugging this estimate into Eq. (1.7) with the 
limit kernel of Eq. ( ^.25 ), we first get 



(1-A) 



(3.28) 



c\ can be determined by mat ching the coefficients of 
t he le ading terms in Eq. ( |l.5| ) using the kernel of Eq. 
( |3.25 ). After a straightforward calculation, one gets in 
the d — > oo limit, 



cx = 2{l-X)Ix 

r-l 



/a = / [u{l - u)]-^-^'^ [u^ + {l-u 



which leads to 



l + \~2^-'^I 



dr-l 



(3.29) 

1^ - l] du (3.30) 



(3.31) 



We thus find a nonperturbative (exponentially small) cor- 
rection to T in the large d and large D limit, consistent 
with the result obtained above for d > 1 and large D. 
Note that Eq. (3.29) is also consistent with the exact 
result that r ^ 1 as Z? ^ oo for finite d > 1, a result 
that we obtain by setting A = (as I\ diverges). 



D. Summary of the results 



C. Large d and D 

We now present a nonperturbative calculation in the 
limit of large d and D, keeping the ratio X = d/D fixed. 
In this limit, the kernel can be written, 

(a;* + yif = 2'^{xy)i{l + 0{d/D^)) (3.25) 

and surprisingly transforms into the well-studied "prod- 
uct" kernel [^| jT^-|2^ , p^ . Assuming scaling (a still contro- 
versial subjcct|21|), one ca n ea sily show that t — 1 + A = 
1 + d/D ^ (see also Eq. (Q and (|oJ) and the discus- 
sion below them, as it corresponds to fi — X/2 > 0). 

We can show that including higher order corrections in 
power oi 1/D does not change the value of r such that 
the correction to r = 1 -I- A is certainly nonperturbative. 
Consider the expansion of the kernel: 



We have shown that when D ^ oo, r — > 1 for d > 1, 
whereas r ^ 2 — 2^~'^ < 1 for d < 1. We were able to 
derive an 0(1/1?^) perturbative expansion in d < 1, and 
we convinced ourselves that the leading corrective term 
in d > 1 was of order 1/D, by giving an estimate of this 
correction. In d = 1 both approaches break down and 
the large D corrections to Too — 1 are probably nonper- 
turbative. 

When d 0, T goes to zero, and we gave a first order 
perturbative expression in d, for any D. For D — 1, we 
also found the explicit coefficient in d^. 

Eventually, we showed that for a fixed homogeneity 
A = d/D, T tends exponentially to 1 -I- A at large d. In 
the following section we present a new general numeri- 
cal method to compute r and we confirm our analytical 
result by performing the first extensive study of the func- 
tion T(d, D). 
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IV. VARIATIONAL APPROACH 

In this section, we present a practical way of obtaining 
good approximate values for r, without explicitly solv- 
i ng S moluchowski's equation. Once again, we rely on Eq. 
(1.8), which holds for the exact scaling function (solution 
of Eq. (1.5)), for any a > r — 1. This equation is general, 
and does not depend on the specific kernel we study in 
this article. As a consequence, the methods we develop 
are general and do apply to any homogeneous kernel. We 
emphasize the fact that this method does not intend to 
approach the whole scaling function, but sets the focus 
on the computation o f r ( in fact, numerically solving the 
scaling equation Eq. (1.5) for the scaling function seems 
to be at least as difficult as directly solving the time- 
dependent equation [p9|). 



A. Principles of the method 

The simplest way of approximating t is to evaluate the 
"average" in Eq. ( |2.2| ) using a reasonable trial weight 
function g{x,y) instead of the unknown exact one. As a 
simple start, we will expose a crude, but straightforward 
algorithm, that illustrates the basic idea. Then we will 
develop the variational method itself, which is not much 
more intricate, but much more effective. 

A one parameter choice for a trial weight function is ob- 
tained by replacing in the above expression of g{x, y) the 
exact f{x) by frix) = x~'^ exp(— x) which has the correct 
leading asymptotics for small x (by definition of r) and 
decays exponentially at large x, although not with the 
exact asymptotics x~^e~^ (A < 1) Still, this func- 

tional form is known to be a good approximation of the 
actual f{x) obtained in simulations pl|], and is even the 
exact solution, but for a multiplicative constant, for the 
constant kernel (r = 0) and in the d = D = 1 case, which 
belongs to the special class A = 1 |^ . The first idea that 
comes to mind is j ust to determine r self-consistently 
such that Eq. (2.2) holds for /,-, with a specific choice 
of a, for instance a = A. This is readily done, by an 
iterative method: starting from an initial tq, verifying 
previously obtained exact bounds, we construct the se- 
quence. 



Tn+l = (1 - £) + e (2 - (1 - a)RMrJ) 



(4.1) 



with 

i?a(0) -2- 



IIo °° x"(l){x)y^(j){y)dxdy 



J /+°° <j,ix)<j,{y)K{x, y) [{x + y)» - x» - 



which converges, with a proper choice of 1 > e > 0, t o 
a fixed point corresponding to an fr verifying Eq. (2.2). 
The numerical evaluation of Rij) can be achieved with 
utter celerity and arbitrary precision, since it reduces to 



the calculation of one-dimensional integrals, and of a few 
values of the gamma function, thanks to a very conve- 
nient transformation (see appendix ^). We notice that 
it is unnecessary to include any multiplicativ e co nstant 
into fr, since it would just cancel out in Eq. ( ^.2| ). 

Of course, this algorithm should yield different values 
of T for different choices of a, except in the special case 
when the exact solution is of the form /,-. This corre- 
sponds to d — Q, D — oo and d = 1,13 = 1, and this 
method converges by construction, to the exact value 
of r, but for the round-off errors. In the generic case, 
the variation can be non negligible (in d = 2,D = 4, 
T « 1.371 for a = d/D, while r « 1.398 for a = 0.403)) 
and the fixed point r may even violate exact bounds. For 
instance, in the case d = 1, 13 = 3 with a = d/ D we get 
r = 0.9894 whereas we know that r > 1. The variation 
with a makes the method unreliable. In d = 2, Z3 = 4, it 
gives T « 1.385 ± 0.015, compared to r « 1.434 ± 0.004 
with the variational approximation, that we now intro- 
duce, which, starting from the same basic idea, proves to 
be much more effective. 

Variational approximation - A much better and hardly 
more intricate method is to choose a reasonable sample 
of values of a, and minimize an error function measur- 
ing the violation of the corresponding Eq. (2.2). This 
method can be systematically improved by allowing for 
n free 'fitting' parameters (including r itself) in the trial 
weight g{x,y). In the following we will proceed by re- 
placing the exact / by a variational function of the form. 



/^(x,ro,ri, ..Tn,ci..Ci) 



Y^c,x-^^e-^ (4.3) 



and we will minimize the error function, 

= E - 2 + (1 - a,)^o, (/.))' 



(4.4) 



to get a variational approximation t.„ = tq of r. Brute 
force should no t be used in the evaluation of once 
again, Eq. (1.1) makes it possible to drastically reduce 
the computation time, and to perform the evaluation of 
with an excellent precision. 

Of course, the values of the exponents in /„ should not 
be blindly chosen, van Dongen and Ernst ||l^ showed 
that the subleading term in the small x asymptotic ex- 
pansion of / is 



2,l + A-2r^ if T > 1 + A-/X1, 

cx \ x^i-^, if T < 1 + A - /xi 
x~'^ In a;, if r = 1 + A — /^i 



(4.5) 



(4.2) with K(x, y) — x^ oc y'^^x^ when x ^ oo, whereas the 



exact asymptotic at large x is oc x e ^. Therefore, a 
good three-parameters class of trial functions should be: 



/l,(x,To,Cl,C2) 



1 



Cl 



C2_ 



(4.6) 
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Ti being either 2ro — 1 — A (if tq > 1 + A — /zi), or tq — /ii 
(if To < 1 + A — ^i). The small x leading term in fy is tq 
provided that tq > X. The approximate value r„ is the 
value of To at the minimum. 

By construction, this method reproduces the exact re- 
sults for the constant kernel and d = 1, D ^ 1, since 
the exact scaling function is contained in those cases in 
the class of variational function we chose. In general, 
this method is inadequate to approach / itself, and is 
just designed to compute r, in the same way as the vari- 
ational approach in quantum mechanics is designed to 
obtain eigenvalues but, in principle, not eigenfunctions. 



B. Implementation 

With a small number n of variational parameters, we 
choose to perform the minimization with the downhill 
simplex method described in |Q (steepest descent, con- 
jugate gradient or other methods could also be used, with 
the drawback that these methods require extra evalua- 
tions of to compute its gradient). This method starts 
from a n-dimensional simplex, i.e. n + 1 points in the n- 
dimensional parameter space, and performs a sequence of 
geometric deformations until it contracts to a local mini- 
mum of the function. It is not the fastest algorithm, but 
it easily converges, and in our case where the computa- 
tional burden is low we do not need more sophisticated 
devices. 

As in any optimization problem, the initial condition 
is a crucial parameter, but here there is the additional 
complication that the smallest moment amin used in the 
computation of x^, should be bigger than r — 1, and 
bigger than tq — 1 at any step of the algorithm. What 
information on the value of r we may a priori gather (ex- 
act bounds, perturbation expansion), should guide our 
choice. Anyway, we do know that t < 1 -I- A: starting 
with an initial tq smaller than 1 + A and amin > A should 
avoid any trouble. As we get a first approximation of r 
we will be able to decrease the value of amin and make 
it closer to Ti, — 1, while refining the initial conditions. 
A few Monte-Carlo minimization steps can also be used 
to find a proper initial condition (but we scarcely needed 
this functionality in this work). 

Why should we choose as small an a^in as possible 
? The answer is that small moments probe the small 
X divergence of f{x), which is precisely what we are in- 
terested in. However, we also need some intermediate 
and higher moments to probe the intermediate x and the 
large x decay to stabilize consistent values of ci and C2. 
There should be at least as many moments as variational 
parameter, otherwise there would be an infinite number 
of minima. Too many moments would cause excessive 
numerical round-off errors in the computation of x^- 

We tested round-off errors by computing t^, for the 
exactly solvable model Kl for which f{x) oc x~^^^e~''^, 
since, were we endowed with infinite numerical precision. 



our algorithm would yield the exact result in this case, as 
said before, whatever the ai may be, provided that they 
all are bigger than 1/2 = t — 1. 

With the three parameter function introduced above, 
and moments 0.55, 0.667, 0.783, 0.9 and 2, we find 
T = 1.49997 ± 4 X 10^*^ (x^ = 1.94 x 10^^), the un- 
certainty being due to variations with different choices 
for the initial values of the parameters and the tolerance 
on the size of the simplex (the minimization algorithm 
stop criterion). The round-off errors increase with the 
number of moments and the number of variational pa- 
rameters. The error is much bigger on ci and C2, we find 
ci = 0.11 ± 0.1 and C2 = -0.12 ± 0.1, instead of strictly 
0. This means that the sensitivity on ci and C2 is small 
in the vicinity of the minimum, and this method is not 
the right one to determine the scaling function (a nega- 
tive C2 is unphysical here), but it just was not devised for 
this purpose: we just meant to compute r, and for this 
quantity the accuracy is excellent. 



C. Numerical results 

We used this method to determine approximations of 
T for the kernel {x^^^ +y^^^)'^ . We compared our results 
to numerical values obtained for d < 1, D = 1 by Kriv- 
itsky [pl| , and to our perturbative and nonperturbative 
expansions. 

All values were obtained from the three-parameter 
variational functions introduced earlier in this text. We 
used 8 moments, 6 in the interval [amin, 0.9], plus a = 2 
and a = 3. amin was adjusted to be as close to Ti, — 1 as 
possible. The computation time was from 1 to 10 seconds 
per run on a HP workstation. 2 to 5 runs per points were 
necessary to adjust the parameters. 

We also computed a few points with a different repar- 
tition of moments: 5 in the range [t — 1, d/D], a = 0.9, 2, 
3, as well as with only 2 variational parameters (ci = 0), 
and with 4 variational parameters (the additional expo- 
nent being /ii — t in the case when t > 1 -I- (d — 1)/ D). 
The observed relative variations of t„ were at most of a 
few 10"'^. In all cases, t was found to be consistent with 
exact bounds. 

First, we consider the case D — \. Fig. (|^) shows 
the comparison between variational approximations of 
T obtained with the modus operandi we just exposed, 
values extracted by Krivitsky [pl|| from a numerical so- 
lution of Smoluchowski's equation, and the 0{<P) per- 
turbative expansion. The agreement between the vari- 
ational approximation and Krivitsky's results is excel- 
lent, which confirms the effectiveness and efficiency of 
the method: the ratio computation time (a few sec- 
onds)/accuracy is impressive. Actually, the variational 
approximation looks smoother than Krivitsky's curve, 
which has two visible accidents (small cusps) near d = 1 
and d = 0.4, and the variational approximation is fully 
consistent with the exact 0{(P) expansion at small d to 
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which it clearly tends asymptotically, whereas Krivitsky's 
result tends to remain parallel to the perturbative curve, 
though close to it. Its good agreement with our infinite 
time results assesses the fact that Krivitsky's solution 
actually reached the scaling regime, which, as said in 
section |, was not obvious a priori. We conclude that in 
this regime, the variational approximation recovers and 
confirms the results obtained by numerical integration of 
Smoluchowski's equation. 



1.6 
1.4 
1.2 



1 numerical (Krivitsky) 
variational 

perturbative 2nd order 




tent with a nonperturbative decay in l/D (see below). 
For d > 1 the large D decay is slower as analytically 
predicted (we found a 1/D perturbative correction, see 
below) . For d>2 the curves shape qualitatively changes 
and an inflexion point appears. 

Iso-A lines exponentially saturate to 1 + A at large £), 
as analytically established before. Fig. ^ shows the com- 
parison between the variational approxima tion and the 
nonperturbative large d expansion of Eq. (3.31) in two 
cases, A = 1/2 and A = 2/3. The agreement is once again 
excellent at large d. 

In d = 1,D = 2, Song and Poland H found r = 
1.123 ± 0.016 (using their first result), which compares 
well with our t = 1.109. In d = 2, D = 3, we find 
T = 1.528 which, unlike their result (1.243), is perfectly 
consistent with the exact bounds 1.4349 < r < 1.585. In 
d = 2, D = A, we find r = 1.434, which is in fair agree- 
ment with the perturbative large D estimate r = 1.462 of 
section III. In fact, as shown on Fig. ^, the perturbative 
estimate is indeed a good approximation of r in d = 2 
for D < 6, and the c>c 1/D decay is confirmed by the 
variational results. The cusp on the variational curve is 
confirmed by the existence of an inflexion point on c? > 2 
curves, as mentioned above. In d = 1, a nonperturbative 
exponential large D decay to Tqo = 1, is confirmed by 
Fig. I We roughly find r - 1 oc e-^-^^^. 



FIG. 1. In D = 1, the comparison between the results ob- 
tained in jil] by Krivitsky, the variational approximation with 
3 parameters and 8 moments, and the O(d^) perturbative ex- 
pansion of r, illustrates the efficiency of the variational ap- 
proximation. Indeed, the agreement between the numerical 
solution of Smoluchowski's equation j2l] and the variational 
approximation is excellent. The variational approximation is 
even in closer agreement with the small d perturbative ex- 
pansion than Krivitsky's result, and although both methods 
recover the exact result t — 3/2 for D = 1, Krivitsky's curve 
seems to have an accident in the vicinity of _D = 1, whereas 
the variational result is smooth. 



Once the effectiveness of the method was established, 
we were able to carry out the first systematic study of 
T(d, D), and to control its validity thanks to the analyt- 
ical results obtained in sections 



Jand^. 

We show on Fig. |, the function r(d, D) (0.25 < d < 3, 
d < D < 8.) plotted in a (r, D) diagram. Two kinds 
of curves are shown. Solid lines represent some iso-d 
lines, i.e the function t(D) for a fixed value of d, whereas 
dashed lines are iso-A (A = d/D) lines. The reliability 
of the approximation is assessed by the comparison with 
analytical results. As established in section [II iso-d lines 
tend to r = 2 — 2^^^* (stars on the right axis of Fig. |^) 
if d < 1, and to 1 if d > 1. As expected, the critical 
D above which r becomes smaller than 1 tends to infin- 
ity when d —^ 1~ , entailing the breakdown of the large 
D perturbative expansion in D > 1. The d = 1 iso-d 
line seems to tend exponentially to 1, which is consis- 
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0.0 



1.0 2.0 



3.0 



4.0 

D 
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FIG. 2. The exponent r was computed by the variational 
method for various values of d and D. We show here some 
iso-d(sohd lines) and iso-A (dashed) (A — d/D) hues. The 
iso-d hnes tend to r = 2 — 2^^^* (stars on the right axis) if 
d < 1, and to 1 if d > 1. The critical D above which r 
becomes smaller than 1 tends to infinity when d ^ 1~, en- 
taihng the breakdown of the large D perturbative expansion 
in D > 1. The d = 1 iso-d line seems to tend exponentially 
to 1, while for d > 1 the relaxation to 1 is slower. An inflex- 
ion point appears above d « 2. The iso-A lines exponentially 
saturate to 1 4- A at large D. 
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► — • variational 
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).0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 

d 

FIG. 3. Iso-A curves computed by the variational method 
(sohd lines), as a function of d, for A = 1/2 and A — 2/3. 
As analytically established, r tends to 1 + A at large D. The 
agreement is good at large d with the nonperturbative expan- 
sion (dashed lines). 



Eventually, we show on Fig. || (for d = 0.25), that 
the variational result is also in good agreement with the 
large D second order perturbative expansion in d < 1 

(oc 1/7^2)^ 

As this section draws to a close, we shall say that this 
variational method, although very simple, seems to be 
very well adapted to the determination of the exponent 
T, as it is fast and, at least in the case we studied in 
this article, very accurate. It made it possible to ac- 
quire for the first time quantitative knowledge of r in 
the whole parameter space of the Kf^ kernel, the most 
studied and the prototype of the notorious class II ker- 
nels. The method is general and could help shedding 
some light on the whole class of kernels, thus increasing 
the practical use of Smoluchowski's approach to under- 
stand aggregation phenomena. This point is worth an 



example. This is precisely what is dealt with in section 
0- 
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FIG. 4. In d = 2, the exponents computed by the varia- 
tional approximation are in good agreement with the pertur- 
bative large D estimate r = 1 -j- 1.849/D. From data, the 
actual asymptotic correction seems to be closer to 1.82/D. 
The cusp on the variational curve corresponds to the change 
of behavior with the occurrence of an inflexion point for above 
d = 2. 
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FIG. 5. For d = 1, the exponents computed by the vari- 
ational approximation displays a much faster decay to their 
D = oo limit {too = 1), than for d > 1. Indeed, as shown 
on this figure, the decay seems to be exponential in D, with 
roughly r — 1 oc e~^'^^^ , a nonperturbative behavior to be 
related to the break-down of the large D perturbative ap- 
proaches for d = 1. 
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FIG. 6. In d = 0.25, the exponents computed by the varia- 
tional approximation are in good agreement with the pertur- 



bative large D estimate t = 2 — 2 + 
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V. APPLICATION IN TWO-DIMENSIONAL 
DECAYING TURBULENCE 

In this section, we would like to illustrate the results 
obtained in this article by presenting an original appli- 
cation outside the field of massive particle aggregation, 
namely the dynamics of vortices in two-dimensional de- 
caying turbulence. 

Recently, a statistical numerical model was introduced 
pO| , pl| which describes the dynamics and the merger of 
vortices with the assumption that the typical core vortic- 



ity oj and the total energy E ^ J 'i 



conserved {Ri is the radius of the i-th vortex) through- 
out the merging processes. This model reproduces the 
main features observed in direct numerical simulations 
(see I|3^j3l|] for details). For instance, after noting that 
a distribution of vortex radii satisfying P{R) ^ R^^ is 
equivalent to a Gaussian energy spectrum E{k) ~ k^^^ 
pl[ , the simulation of this model was able to repro- 
duce the fact that starting from a Batchelor spectrum 
E{k) ~ {P = 3), the system evolves systematically 
to a steeper spectrum E{k) ~ k~'^ with 7 = 6 — /3 in the 
range 7 « 3 '--^ 5 |3l] . 

Now, one expects that the collision kernel between two 
vortices is somewhat intermediate between the ballistic 
hard-disk form a ^ + R2) and the totally un- 
correlated form a ^ (i?i +i?2)^ (where the probability of 
colliding is proportional to the probability that tw o ra n- 
domly placed vortices overlap, see also below Eq. (1.1)). 
Thus, one can describe approximately the decay of vor- 
tices due to mergers by means of Eq. (^^) with 1 < d < 2 
and = 4, as two colliding vortices merge into a new 
one with R = (i?f -f- i?|)^^'' in order to conserve energy 
and core vorticity. One thus expects a power law radius 



distribution P{R) ~ R'^ , with /3 = D{t - 1) + 1 and 
T given by our model. We find values of 7 ranging from 
7 w 3.26 for d = 2 (taking r = 1.434) to 7 « 4.95 (taking 
T = 1.012) for d = 1, in good qualitative agreement with 
observed exponents. As also found in direct simulations, 
the actual exponent (and here the value of the effective 
correct d) could depend on the actual initial conditions 
{uj, area occupied by the vortices ^ enstrophy). Note 
that the Batchelor limit case 7 = 3, is obtained when 
taking the naive strict upper bound t ^ I + d/ D with 
d = 2 and £) = 4. 



CONCLUSION 



In this article, we tackled the notoriously difficult 
problem of nontrivial polydispersity exponents in Smolu- 
chowski's approach to aggregation from an original angle. 
We chose to directly start from the scaling (infinite time 
limit) equation, and we did not focus on the determina- 
tion of the whole scaling function, which is the object 
of solving Smoluchowski's equation, to concentrate on r 
itself, which actually mainly depends on global (integral) 
equations. We think, and illustrated this point on the 
example of a simplified model of two-dimensional turbu- 
lence, that in some cases, the only knowledge of r would 
still be a good step towards the understanding of the phe- 
nomenon. The choices we made were fruitful and gave 
birth to new analytical and numerical results. 

From an analytical viewpoint, we were able to use in- 
tegral equations to find some exact bounds for r, and, in 
the specific case of K'^ = (x^/^ -I- y^^^)'', we obtained 
some perturbative and nonperturbative expansions of r, 
without explicitly computing the corresponding expan- 
sions for the whole scaling function. 

From a numerical viewpoint, we devised a variational 
approximation scheme, that recovers by construction 
known exact results, and can be used as a tool for ex- 
tensive determination of r, since it is both very econom- 
ical and accurate. In addition, it is likely that the scal- 
ing function obtained in the variational approach is in 
many cases qualitatively, if not quantitatively, right. To 
illustrate its effectiveness, we performed a comprehen- 
sive study of r for a wide range of the parameters (d, D) 
of the kernel Kf^. This is a noticeable advance, since 
very little quantitative knowledge was available for this 
kernel, although it was the prototype kernels with a non- 
trivial t, and the object of much attention in the past 

©lljlolH. 
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APPENDIX A: A USEFUL FORMULA 



J 



where T is the gamma function, and, 

K{l,u) [(1 + w)" - 1 - ■ 



x-''^y-''^e-''-yK{x, y) [{x + y)" - - y"] dxdy = 
r(2 + A + a - Ti - T2) [X{Ti,a,Ti + T2) + X{T2,a,Ti + T2)] , 

I 



(1.1) 



X{t,a,q) = 



u*(l + u)2+^+"-9 



Uu (1.2) 



It would be very awkward and inefficient to use 2- 
dimensional numerical integration (especially here, as the 
integrand is singular at the origin). A startlingly eco- 
nomical way of computing the gamma function is due to 
Lanczos and is described in p3] (it is not much slower 



To demonstrate this formula is straightforward: just than the built-in exponential function...), 
make the change of variable x = uv, y = v, and use the 
definition of the F function: 



Fix) 



(1.3) 



APPENDIX II: THE 0{D^) TERM IN D=l 



From a numerical viewpoint this formula makes it We derive the 0{d^) correction to t ^ 2d fo r D — 1, 
possible to implement very rapid and accurate code by c omputing the d^ order of respectively Eqs. (1.7) and 
for the variational approximations we developed before. (1.8) with a — d, to get, 



4 - 2a2 c -I- 4 



fi{x) In xdx + A / f2{x)dx 



-t-00 



-02 = - 

2 



+ 00 



+c>o 



e ^(Inx) dx + fi{x) In xdx -\- / /2(x)dx, 



(2.1) 
(2.2) 



where t = 2d + 026?^ + 0{d'^), and. 



c = 4 



+ 00 



e-^(ln(x))^da; + 4 



+00 



-x-y 



ln(x -|- y) In 



xy 
x + y 



dxdy 



+00 



fi (x) dx 



e y \ny dy 



fi (x) dx 



(2.3) 



r 



c can be computed since J fi is known from the first the small d expansion of the scaling function. With 
order calculation. After some elementary transforma- u{x) = e^/i(x), the latter equation is 
tions, we find that c — 4 f (\nx)^dx — 



Combining Eq. (^ and (gj)), we find 4 + 2a2 
c — 4 / e^^{lnx)'^dx, hence eventually 

^2 



xu" + (1 — x)u' ^ ^ 



+00 



a2 = y-4. 



(2.4) 



D 



\nx — 2Jd- 



(3.1) 



APPENDIX III: THE LINEARIZED SCALING 
FUNCTION 



With v{x) = u{x)/{x — l), this equation reduces to a first 
order differential equation for v' , and we find. 



We find the solut ion of the second order differential 
equation Eq. ( 3.10 ) for the linear coefficient fi{x) in 



fi{x) = Co uo{x)e ^ + ci{x - 1) - 2Jd - -^(1 + Inx) 



dyi 



2/1 (yi 



-2/ dy2y2 ey^y2~l) dy^ 

1 ~ ^) Jo Jo y^' + y,' 



(3.2) 
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and 



APPENDIX IV: PERTURBATIVE ESTIMATE 



uo{x) = e"^ - {x - l)Vp 



y 



-dy 



(3.3) 



( "Vp" means "principal value" ) . 

In fact, the triple integral can be transformed into a 
simple integral involving special functions. For our pur- 
pose, we only need to know that this integral goes to zero 
when X —> 0, which is easily seen. 



For d > 1, T = 1 + e{D) where e ^ when D 
We make the ansatz: 

fix)^foo{x) + 



(4.1) 



and plug it into Eq. (1.7) to obtain, 1 — e = 2^ 
cT{jj — e), which means that, when D — > oo, c 



(1 - 2^-'^){d/D ~ e). Then we make use of Eqs. ( [li| ) 
and (1.8) to obtain. 



2(1 - - e) = ^""'"^ / / e-"-^(a;* + y^f[xT^ + ~ {x + y)*]dxdy 



D' 

+ 2^-'^c r(l + 2d/D - e) [X{Q, d/D, 1 + e) + A:'(l + e, d/D, 1 + e)] 
+ 2c^T{2d/D - 2e)X{l + e, d/D, 2 + 2e) (4.2) 



The next step is to write down the limit of this equa- 
tion when D — > oo. We know that r(a;) '^x^o ^/x, and 
a change of variable v — u'^/'^^^ in the integral factors X 
shows that X{1 + e, d/D, 1 -h e) ~ ^"(1 -I- e, d/D, 2 + e) - 
{d/D - £)-! (1 -I- v^-Ydv. We obtain: 



2 = 2 



2-d 



d/D -S Jq 



{l + i—fdv 



{d/D 



{l + v^Ydv 



(4.3) 



HL is the limit of De. Taking into account the value of c, 
we finally get, 



-tl-d 



{1 + v — fdv = J(k, d) 



- = i + ^ + o(4) 



D 



'D^ 



(4.4) 
(4.5) 



The equation Eq. (4.4) has a unique solution < k < c? 
since the integral J(k, d) is a decreasing function of k, and 
J(0, d) = 2'^> 2/(21-'' _^ 1) > 1 ^ fj^j (for ^ > i). 
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